International Journal of Nanomedicine 



Q Open 



Access Full Text Article 



Dovepress 

open access to scientific and medical research 

ORIGINAL RESEARCH 



Heuristic modeling of macromolecule release 
from PLGA microspheres 



Jakub Szl<jk' 
Adam Paclawski 1 
Raymond Lau 2 
Renata Jachowicz 1 
Aleksander Mendyk' 

'Department of Pharmaceutical 
Technology and Biopharmaceutics, 
Jagiellonian University Medical 
College, Krakow, Poland; 2 School of 
Chemical and Biomedical Engineering 
Nanyang Technological University 
(NTU), Singapore 



Video abstract 



1 




Point your SmarlPhone at the code above. If you have a 
QR code reader the video abstract will appear. Or use: 
http://dvpr.es/H8U6aG 



Correspondence: Aleksander Mendyk 
Department of Pharmaceutical 
Technology and Biopharmaceutics, 
Jagiellonian University Medical College, 
Medyczna 9 St, 30-688 Krakow, Poland 
Tel +48 12 620 56 04 
Fax +48 12 620 56 18 
Email mfmendyk@cyf-kr.edu.pl 



This article was published in the following Dove Press journal: 
International Journal of Nanomedicine 
29 November 2013 

Number of times this article has been viewed 



Abstract: Dissolution of protein macromolecules from poly(lactic-co-glycolic acid) (PLGA) 
particles is a complex process and still not fully understood. As such, there are difficulties in 
obtaining a predictive model that could be of fundamental significance in design, development, 
and optimization for medical applications and toxicity evaluation of PLGA-based multiparticulate 
dosage form. In the present study, two models with comparable goodness of fit were proposed for 
the prediction of the macromolecule dissolution profile from PLGA micro- and nanoparticles. 
In both cases, heuristic techniques, such as artificial neural networks (ANNs), feature selec- 
tion, and genetic programming were employed. Feature selection provided by fscaret package 
and sensitivity analysis performed by ANNs reduced the original input vector from a total of 
300 input variables to 2 1 , 17, 16, and eleven; to achieve a better insight into generalization error, 
two cut-off points for every method was proposed. The best ANNs model results were obtained 
by monotone multi-layer perceptron neural network (MON-MLP) networks with a root-mean- 
square error (RMSE) of 15.4, and the input vector consisted of eleven inputs. The complicated 
classical equation derived from a database consisting of 17 inputs was able to yield a better 
generalization error (RMSE) of 14.3. The equation was characterized by four parameters, thus 
feasible (applicable) to standard nonlinear regression techniques. Heuristic modeling led to the 
ANN model describing macromolecules release profiles from PLGA microspheres with good 
predictive efficiency. Moreover genetic programming technique resulted in classical equation 
with comparable predictability to the ANN model. 

Keywords: poly(lactic-co-glycolic acid) (PLGA) microparticles, genetic programming, feature 
selection, artificial neural networks, molecular descriptors 

Introduction 

Poly(lactic-co-glycolic acid) (PLGA) microparticles play a prominent role in many 
drug formulations. PLGA is flexible and has a satisfactory safety profile. It also has 
the ability to modify the drug dissolution profile for controlled-release formulations. 1 
Moreover, the drug protective properties of PLGA were found to be suitable for unstable 
active pharmaceutical ingredients (APIs). 2 Therefore, PLGA has wide applications 
as an excipient. 

PLGA can be formulated as films, discs, microcapsules, and nano- or micro- 
spheres. 3 " 6 The current study is focused on the formulation of PLGA nano- and micro- 
spheres. One of the most studied properties of PLGA particles over the past years is the 
API release mechanism. There are often contradictory results reported in the literature, 
which demonstrate the high level of complexity of PLGA-based dosage forms. The drug 
release from the PLGA matrix is mainly governed by two mechanisms: diffusion and 
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degradation/erosion. 7 The work of D'Souza and DeLuca 8 and 
Mollo and Corrigan 9 showed that the drug release profile can 
be divided into two stages. Initially, the release is considered 
to be diffusion-controlled, corresponding to the low amount 
of released drug. Afterwards, the drug release is controlled 
primarily by the degradation and erosion of the PLGA matrix. 
There are many factors influencing the diffusion and degrada- 
tion rate of PLGA. For example, pore diameters, matrix-API 
interactions, API- API interactions, and formulation composi- 
tion; 10 " 13 however, there is, to date, no versatile quantitative 
model describing the relationships among these factors in a 
consistent, mathematical manner. Recently, Fredenberg et al 
proposed four so-called "true" drug release mechanisms 
from PLGA matrices: 1 ) diffusion through water- filled pores; 
2) diffusion through the polymer; 3) osmotic pumping; and 
4) polymer erosion. The release profile of large and hydro- 
philic molecules, such as proteins and peptides, is mostly 
limited by diffusion through water-filled pores. 7 

Over the past 2 decades, there have been several publi- 
cations covering the model development of the drug release 
profile from PLGA microspheres. Among the most recog- 
nized are the works of Zygourakis and Markenscoff 14 and 
Gopferich, 15 in which Monte Carlo and cellular automata 
microscopic models were introduced. The authors used a 
two-dimensional model of the polymer matrix based on the 
physicochemical characterization to predict erosion and 
therefore the release rate of small molecules. Later, Siepmann 
et al implemented a model coupling a Monte Carlo simula- 
tor with sets of partial differential equations. 16 The model 
described the chemical reactions and physical mass transport 
processes involved in erosion-controlled drug release from 
PLGA beads. Nevertheless, the predictability of classical 
models (empirical/semiempirical or mechanistic models) 
has not been thoroughly assessed in many studies and the 
versatility remains questionable. 14-16 A recent study by Barat 
et al 17 applied Gopferich 's theory 15 on multi-agent systems. 
Good agreement between modeling and experimental results 
was obtained for polymer erosion coupled with a partial dif- 
ferential equation model of lysozyme release profiles from 
PLGA spheres; 17 however, to the authors' best knowledge, 
there is no general model able to predict the release rate of the 
various macromolecules from PLGA systems. Therefore, 
the design of the PLGA microsphere formulations loaded 
with macromolecules depends heavily on the trial-and-error 
approach. 

The obj ective of this work was to demonstrate the cooperative 
use of feature selection, artificial neural networks (ANNs), and 
genetic programming (GP) to create a simple predictive model 



of the macromolecule dissolution rate from the PLGA dosage 
form, based on the gathered literature data. 

Materials and methods 

Data set 

Data were collected from the literature. After careful screening 
of about 200 publications, release rates of 68 PLGA formula- 
tions from 18 publications were selected for data extraction 
(Supplementary material Table SI). Selected PLGA formula- 
tions were prepared using four methods: 1) a double-emulsion 
water-in-oil-in- water solvent extraction process; 2) a solid-in-oil- 
in-water emulsion method; 3) an oil-in-oil solvent evaporation 
technique; and 4) a spray-drying process. Formulations included 
14 model substances: bovine serum albumin; recombinant 
human erythropoietin; recombinant human epidermal growth 
factor; lysozyme; recombinant human growth hormone; hen 
ovalbumin; human serum albumin; beta-amyloid; insulin; 
recombinant human erythropoietin coupled with human serum 
albumin; L-asparaginase; bovine insulin; alpha- 1 antitrypsin; 
and chymotrypsin. The release profiles presented in the literature 
were manually digitized, taking mean values of the observations 
if possible. Standard deviations and/or errors were not included 
due to the manual character of the digitization and lack of this 
information in some sources. Overall there were 745 data 
records with 320 variables (Supplementary material Table S2). 
The independent parameters contained 319 inputs covering the 
formulation characteristics (PLGA inherent viscosity, PLGA 
molecular weight, lactide-to-glycolide ratio, inner and outer 
phase Polyvinyl alcohol (PVA) concentration, PVA molecular 
weight, inner phase volume, encapsulation rate, mean particle 
size, and PLGA concentration); the experimental conditions 
(dissolution pH, number of dissolution additives, dissolution 
additive concentration and production method, and dissolution 
time); and the molecular descriptors of the macromolecules 
and excipients. The molecular descriptors were computed using 
Marvin cxcalc plugin, UK (version 5.11; ChemAxon, Budapest, 
Hungary) for the drug substance as well as the excipients. 18 The 
amount of the drug substance released (Q) was the only depen- 
dent variable. Prior to feature selection, the initially obtained 
319 inputs were further reduced to 300 inputs by removing 
all the null and missed input during the data set preprocessing 
procedure. The data set was then processed in several ways: 

• Noise addition to prevent models from over-fitting. The 
noised data records were produced numerically with +5% 
amplitude for each variable value and two times more 
records number (noted in the text as "rand"), 

• Data set split according to the tenfold cross-validation 
scheme, with the aim of excluding all the data belonging to 
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the particular formulation, thus simulating the real application 
of the model forced to predict the behavior of the unknown 
formulation (noted in the text and tables as "lOcv"), 
• Linear scaling using either the output range of <0.2, 0.8> 
or <— 0.8, 0.8>, in order to match nonlinear activation 
functions domains were used for ANNs only (noted in 
text and tables as "scale"). 

Feature selection 

It is considered state of the art to implement the feature 
selection technique prior to the actual modeling stage when 
dealing with multidimensional data sets. Two variable selec- 
tion techniques were used to reduce redundant or irrelevant 
information gathered during data set preparation: based on 
the methodology proposed by Zurada et al 19 with further 
modification by Mendyk and Jachowicz, 20 which uses ANNs 
trained with the back propagation algorithm as the model- 
ing tools; and feature ranking created by the fscaret pack- 
age of the R environment (The R Foundation for Statistical 
Computing, Vienna, Austria). 21 - 22 

The main parameters of the applied techniques are listed 
in Table 1 . The resulting feature rankings were then inspected 
on the substantial decrease of variable significance; two cut- 
off points were chosen for each method. 

ANNs 

Based on reduced vectors of inputs, two types of ANN model 
were obtained, namely, multi-layer perceptron artificial 
neural networks (MLP-ANNs) and monotone multi-layer 
perceptron artificial neural networks (MON-MLP) networks. 
The tools used during modeling were the Nets2012 ANN 23 
simulator (developed in-house, Aleksander Mendyk, 
Jagiellonian University, Krakow, Poland), and R environment 
with the monmlp 24 package for MON-MLP networks. 23 ' 24 



Table I Feature selection parameters used in the assay 
according to applied methods for artificial neural networks and 
fscaret 



ANNs 



fscaret 



Feature selection technique applied 

Whole data set. Whole data set. 

Trained and tested over Time-limiting function was set 

170 ANN architectures. to 2 hours for single model 

5M iterations. development. 

Final ensemble based on the "PreprocessData" function was 

goodness of fit criterion: up to off. 

200% of minimum error achieved. The feature ranking was created 

according to the results of 
58 models. 

Abbreviations: M, millions; ANN, artificial neural network. 



MLP-ANNs and MON-MLP were trained on reduced vectors 
after feature selection procedure. 

MLP-ANNs 

MLP-ANNs used the backpropagation training algorithm 
and were tested with activation functions such as linear 
("lin"), logistic ("sigma"), hyperbolic tangent ("tanh"), and 
logarithmic function ("fsr"). The architecture of MLP-ANNs 
comprised one to seven hidden layers. In addition to the 
MLP-ANNs, neuro-fuzzy systems (NFs) were employed by 
adapting the simplest Mamdani type and contained from five 
to 100 nodes in the hidden layer. All models were multiple 
input/single output (MISO) type. 25 Considering a tenfold 
cross-validation, a total of 3,680 models were trained in a 
single run. Each model was trained up to 10,000,000 epochs 
with several stop points. The generalization error was assessed 
in order to find the most optimal training conditions. The stop 
points were selected as 100,000, 200,000, 500,000, 1,000,000, 
1,500,000, 2,000,000, 3,000,000, 5,000,000, and 10,000,000 
epochs. The epoch size was equal to 1 . Other parameters were 
as follows: 

• momentum technique, with the momentum factor 0.3, 

• delta-bar-delta algorithm, with the initial learning factor 
0.65, 

• jog-of-weights technique designed to prevent getting 
stuck in the local minima of the cost function; a simple 
noise addition to the weights was performed when the 
ANN was not improving its efficiency during the 1 00,000 
epochs (the patience criterion). 

The root-mean-square error (RMSE) was used to mea- 
sure the goodness of fit and was calculated according to 
Equation 1; additionally, relative RMSEs (relRMSEs) for 
each dissolution profile were calculated (Equation 2). 



RMSE- 



i 



n 2 

( pred — obs j 



[i] 



where obs., pred. = observed and predicted values respec- 
tively, i = data record number, and n = total number of 
records. 



RMSE 

relRMSE = ^-100%, 

' e -0„„ 



[2] 



where RMSE = RMSE for formulation i and O and 
Qmin = maximum and minimum percentage of released mol- 
ecule, respectively. 
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MON-MLP 

MON-MLP networks were created and trained by the 
"monmlp.fit" function of the monmlp R package. 24 Bootstrap 
aggregation was used to avoid over-fitting of models. All 
models consisted of two hidden layers with three to 32 nodes 
per layer. Two transfer functions for hidden layer modifica- 
tion of hyperbolic tangent ("tansig") and for output layer 
linear function were applied in each model. The prepared 
ensemble system consisted of ten or 20 neural networks. 
Performance of the models was assessed using RMSE. Other 
parameters were as follows: 

• "trials" function, set to 5 , was used to avoid local minima 
in ensemble, 

• Iteration number for learning models was set to 1 00, 500, 
1,000, or 2,000. 

GP and symbolic regression 

Mathematical equations were produced by means of GP and 
the symbolic regression mode available in the rgp package of 
the R environment. 26 Direct mapping of the input and output 
variables was applied. GP was performed on reduced input 
vectors after the feature selection procedure. 

The parameters of the function "symbolicRegression" of 
the rgp package were set as follows: 

• "individualSizeLimit" (the chromosome length) varied 
from 10 to 300 and was the subject of the optimization 
in regard to the RMSE (Equation 1), 

• "populationSize" was set to 10,000, 

• "myfunctionSet," a set of prototype functions, was restricted 
to the simplest arithmetic operators like addition, subtrac- 
tion, multiplication, and division, together with power, 
natural logarithm, square root, and exponent function, 

• the algorithm stop condition makeFitnessStopCondition 
(RMSE) was equal 1.0, 

• "makeTimeStopCondition" was set to 1 hour for the 
indirect and 120 hours for the direct modeling mode. 
The direct mapping approach applied GP on the origi- 
nal data set to create functional relationships between the 
amount of the macromolecule released (Q) and the vector 
of parameters containing the time variable (t), formulation 
characteristics, and molecular descriptors. 

Although GP provides the model with optimized values of its 
parameters, to achieve better generalization, the multivariate opti- 
mization of equation parameters was applied using the optimx 
package (The R Foundation for Statistical Computing). 27 

Parameters used for the "optimx" function were: 

• "par" - initial value for the parameters was set to 0. 1 , 



• "all.methods" was set to "true", which means that all 
available (and suitable) methods were used 

• "follow.on" was set to "false." All optimization methods 
were used separately. 

Fitting was performed in the mode of tenfold cross- 
validation. The latter procedure was done with the same 
training/testing data sets as for MLP and MON-MLP ANNs 
in order to compare the two systems' generalization abilities. 
The selection criterion for the parameters was the model 
goodness of fit expressed as the RMSE obtained on the 
training data set. 

The outline of the workflow, in conjunction with some of 
the results of feature selection, is presented in Figure 1 . 

Hardware environment 

All calculations were performed using 21 professional 
workstations equipped mainly with Xeon central process- 
ing units (Intel Corporation, Santa Clara, CA, USA) with 
a total of 120 cores and at least 16 GB of RAM operating 
on the openSUSE 12.3 (x86_64) operating system (SUSE, 
Nttrnberg, Germany). 

Results and discussion 

Feature selection 

According to the feature ranking created by ANNs and 
fscaret, four reduced input vectors containing 21,17,16, and 
eleven independent variables were selected. The resulting 
input vectors for 17 and eleven independent variables, which 
yielded the best results for GP and ANNs, respectively, are 
shown in Tables 2 and 3. They originate from two different 
feature-ranking methods and gave the lowest RMSE when 
applied to ANNs and GP. Evaluation of the input vector was 
based on the generalization error (RMSE). A substantial 
reduction of RMSE from 22.8% (ANN models trained on 
tenfold cross-validation with all 300 inputs) to 15% — 1 8% was 
obtained for all four reduced input vectors (Tables 4 and 5). 
Further attribute reduction applied to reduced vectors was 
unsuccessful, giving about 2%-3% higher generalization 
errors in each case. 

Reduced input vectors mostly consisted of formulation char- 
acteristics, which is consistent with previous publications stating 
that the method of preparation and parameters describing the 
particles, influence in great extent the drug dissolution rate from 
the PLGA particles. 1,4 ' 7 Sandor et al stated that the release rate of 
proteins is dependent on the protein molecular weight; 28 however, 
distance-based geometry properties coded as hyper- Wiener index 
and Szeged index were not discussed in the literature. 
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Nets2012 
MLP/NFs nets 
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58 models (linear, 
non-linear) 
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17 
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Figure I Schematic representation of the experimental setup. 

Notes: Data acquisition (Phase A). Data preprocessing (Phase B). Modeling (Phase C). 

Abbreviations: ANNs, artificial neural networks; GP, genetic programming; MLP, multi-layer perceptron artificial neural networks; MON-MLP, monotone multi-layer 
perceptron artificial neural networks; NFs, neuro-fuzzy systems; PLGA, poly(lactic-co-glycolic acid). 
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Table 2 Reduction to 17-variable input vector 


V3.ri3.ble 


Description 


Group 


name 






In I 


Hyper-Wiener index 


Macromolecule 


InJ 
\UL 


Szeged index 


UCiLI IULUI i 


h3 


PLGA inherent viscosity (dL/g) 


Formulation 


M 


PLGA molecular weight (Da) 


characteristic 


InS 


PVA inner phase concentration (%) 




In6 


PVA outer phase concentration (%) 




Inl 


Encapsulation rate (%) 




ino 


Mean particle size (|lm) 




In 9 


PLGA-to-plasticizer ratio 




In 1 0 


Dissolution pH 




In 1 1 


Production method: 1) w/o/w; 






2) s/o/w; 3) s/o/o; 4) spray-dried 




in 1 z 


Hyper-Wiener index 


Plasticizer 


in l j 


log u at pn u 


descriptors 


In 14 


log D at pH 1 




In 15 


log D at pH 14 




In 1 6 


Basic pKal 




ln\7 


Time (days) 


Formulation 






characteristic 


Q 


% of macromolecule dissolved 


Output 


Abbreviations: PLGA, poly(lactic-co-glycolic acid); PVA, Poly(vi 


nyl alcohol); s/o/o, 


oil-in-oil solvent; s/o/w, solid-in-oit-in-water; w/o/w, water-in-oil-in-water. 


Table 3 Reduction to eleven-variable input vector 


Variable 


Description 


Group 


name 






Inl 


Szeged index 


Macromolecule 


Inl 


Pi 


descriptors 


In3 


Quaternary structure of 






macromolecule: 1) monomer; 2) dimer 




In4 


Lactide-to-glycolide in polymer ratio 


Formulation 


In5 


PVA inner phase concentration (%) 


characteristics 


In6 


PVA outer phase concentration (%) 




In7 


Encapsulation rate (%) 




In8 


Mean particle size (|Xm) 




In9 


Dissolution pH 




InIO 


Production method: 1 ) w/o/w; 






2) s/o/w; 3) s/o/o; 4) spray-dried 




Inll 


Time (days) 




Q 


% of macromolecule dissolved 


Output 


Abbreviations: PVA, Poly(vinyl alcohol); s/o/o, oil-in-oil solvent; s/o/w, solid-in-oil- 


in-water; w/o/w, water-in-oil-in-water; pi, isoelectric point. 




Table 4 Artificial neural network-MLP results, according to the 


lowest root-mean-square error SE obtained 




Input 


Architecture Iterations Relative 


vector 




RMSE (%) 


2 1 in 


120 SO 30 20 10 tanh scale 3,000,000 


16.9 


I7in 


7_5_3_fsr_rand_scale I0,000,00C 


i 18.0 


I6in 


200 80 40 20 tanh rand 3,000,000 


17.9 


1 lin 


60_20_8_tanh_rand 5,000,000 


16.2 


Abbreviation: MLP, multi-layer perceptron; RMSE, root-mean-square error. 



Table 5 MON-MLP results, according to the lowest root-mean- 
square error SE obtained 



Input 


Architecture 


Ensemble 


Iterations 


Relative 


vector 








RMSE (%) 


21 in 


8_6_trials_S_orig 


20 


1,000 


17.7 


I7in 


1 2_6_trials_S_orig 


20 


2,000 


17.1 


I6in 


1 2_6_trials_S_orig 


10 


2,000 


16.4 


1 lin 


8_6_trials_S_orig 


20 


500 


15.4 



Abbreviations: MON-MLP, monotone multi-layer perceptron artificial neural 
networks; RMSE, root-mean-square error. 



ANNs 

The best result obtained using the ANN modeling tools had 
a relRMSE of 1 5.4% for MON-MLP networks with a vector 
consisting of eleven inputs (Tables 3 and 5). 

Although the generalization error of the ANN models was 
about 15%— 18%, the relRMSE (related to single formulation) 
of the best MON-MLP model range was from 1 .75% to 36.03% 
between formulations, indicating that there is good potential 
for further improvements. To accomplish this task, the data set 
should be enriched with other crucial parameters influencing 
the dissolution, such as the porosity of particles, pore dia- 
meters, chemical descriptors of polymer, and particle size of 
the excipients. Unfortunately, data were not available or scant 
in the selected publications, therefore not useful for systematic 
quantitative analysis. As a consequence, the model has inferior 
performance when compared to classical models, although it 
is important to note that all the predictive model generalization 
errors were a result of the tenfold cross-validation procedure, 
which is not the case in classical approaches. 7-91415 Moreover, 
our models are universal in terms of the drug encapsulated in 
PLGA microspheres, whereas most of the classical models are 
focused on a limited and/or predefined number of drugs. 1 - 4 ■ 5JS 

Symbolic regression 

Obtained using GP classical equation, though complicated 
and derived from a database consisting of 17 inputs, not 
symbolic regression but equation yielded better generaliza- 
tion error (14.3%) and narrower relRMSE, ranging from 
4.76% to 32.83% between formulations, in comparison to 



Table 6 The results of the formulation-to-formulation relRMSE (%) 



Descriptive 
statistics 


relRMSE (%) 




MON-MLP model 


GP model (Eq 3) 


Min 


1.75 


4.76 


Max 


36.03 


32.83 


Mean 


14.39 


13.75 



Abbreviations: GP, genetic programming; MLP, multi-layer perceptron artificial 
neural networks; MON-MLP, monotone multi-layer perceptron artifical neural 
networks; relRMSE, relative root-mean-square error. 
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(Eq. 3) 



MON-MLP 
models 



^ 75 
O 



R 2 =0.9407 




O 














relRMSE =10.2 



25 50 75 

Observed (Q[%]) 



c 
c 
o 

1 

3 

E 
o 




25 50 75 

Observed (Q[%]) 



B 

o 

— 20 



2 =0.4960 










0 ° °°» 


o 


O 








relRMSE 


=11.7 



10 20 

Observed (Q[%]) 



o 
c 

c 
o 

_ra 

3 

E 



O 

lj 10 




10 20 

Observed (Q[%]) 



£=; 75 

o 

■D 50 



R 2 =0.7998 




° 












O / 






O 


relRMSE =5.3 


) 


10 


20 30 40 

Observed (Q[%]) 


R 2 =0.9698 





















relRMSE =5.5 


) 


25 


50 75 



o 
c 

c 
o 
'■*-» 

3 



If) 

CO 

o 
c 
c 
o 
'■*-» 

3 

E 




10 20 30 

Observed (Q[%]) 



S5 75 

o 



50 



relRMSE =6.9 



Observed (Q[%]) 



25 50 

Observed (Q[%]) 



75 

as 
o 

■n 50 



R 2 =0.8364 ^. 






K O < 


ro o 










relRMSE =6.8 


0 25 




50 


75 



to 

o 
c 

c 
o 

re 

3 

E 




Observed (Q[%]) 



25 50 

Observed (Q[%]) 



Figure 2 Comparison of predicted versus observed values for two best models (MON-MLP and classical equation). 
Note: Data from the gathered database (Table S2). 

Abbreviations: GP, genetic programming; MLP, multi-layer perceptron artifical neural networks; MON-MLP, monotone multi-layer perceptron neural networks; relRMSE, 
relative root-mean-square error. 
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the MON-MLP model. The equation was characterized by 
four parameters (Equation 3). Moreover, during the evolution 
process, the GP algorithm further reduced the number of 
necessary variables by eliminating input numbers 2, 4, 9, 10, 
12, 13, 15, and 16. The simplified mathematical model still 
retained all the types of input variables (with one exception 
of the production method), yet its predictive performance was 
comparable to the more complex MON-MLP model (Tables 6 
and Supplementary material Table S3 [full results]). 

Vlnl7 ^Inll+a+e™ 

C2-4C1 

4ln7 ■\n{ln7)-Inl7-[\+ \[m$-C4\ [3] 

InllnS 

In3 + \n(ln6) In11 -^-^ 
1 ' InSInl4 

Inlle 

in which Inl—Inl7 correspond to labels in Table 2 
and C1-C4 = constants. Equation parameters: 
CI =0.1382595046751, C2 =0.1063092562504, C3 
=4.3209748739209, and C4 =0.3010676275288. 

Due to the lack of data, only 68 formulations and 
14 model substances were included. Nonetheless, since our 
model is based on the molecular descriptors, it is able to pre- 
dict the behavior of a completely unknown drug substance. 
An example of such accurate prediction is formulation 54 
(MON-MLP model relRMSE 7.6%, GP relRMSE 10.2%, 
molecule recombinant human erythropoietin coupled with 
human serum albumin [Figure 2A]). 

In general, when considering the relRMSE of formulations, 
almost 60% of the release profiles had a relRMSE less than 
15% when the MON-MLP model was applied. Moreover, the 
equation obtained by GP yielded 64% of profiles with a rel- 
RMSE less than 15%. Formulation-to-formulation comparison 
indicates that ten formulations yielded a relRMSE below 1 0%. 
It is important to note that those formulations consisted of 
medium-to-large molecules with a molecular weight ranging 
from 14.2 kDa (lysozyme) to 66.5 kDa (bovine serum albumin). 
Formulations that have a relRMSE of more than 30% could 
be perceived as outliers, with two formulations consisting of a 
relatively small molecule of bovine insulin (5.8 kDa). 

It has been observed that, in some cases, the GP model 
yielded better predictions than MON-MLP, and vice versa 
(Figure 2B and C). The GP model had over- or underes- 
timated the initial burst of molecule release (Figure 2D 
and E) and the end-point release (Figure 2D), while the 
midpoints of release profiles were in good accordance with 
the observed data. In contrast, for the MON-MLP models, 



even if it yielded worse generalization error, the shapes of 
the predicted profiles were preserved for most formulations 
(Figure 2B-E). 

Based on the exploratory analysis of the modeling 
results, it could be concluded that the GP models are 
inferior in the presence of recombinant human epidermal 
growth factor and human serum albumin, showing much 
higher errors (Supplementary material Tables SI and S2), 
whereas the MON-MLP model responded poorly to the 
presence of beta-amyloid L-asparaginase and bovine insulin 
(Figure 2C). Moreover, if the value of the variable "mean 
particle size" is above 20 |im, the GP model predictions 
are most likely poor (Figure 2B) and in the case of the 
MON-MLP model, the values below 1 |im yielded worse 
predictions. The above presents limitations of our model, 
which will be a subject of its development in the future. 
Due to the heuristic nature of the models, it is difficult to 
explain all the unexpected behaviors of the models, such as 
the erratic profile seen in Figure 2E, where the MON-MLP 
model creates a completely opposite profile at the time of 
20 days. Thus we wanted to present our analysis of the mod- 
els and their advantages (generic for the molecule) versus 
disadvantages (limitations of predictions) in order to allow 
every user to understand their applicability. The predicted 
versus observed release profiles of discussed formulations 
are presented in Figure 2, with the "observed" data derived 
from the literature. 29 " 33 

Conclusion 

In this paper, a new approach was presented for modeling 
the dissolution of macromolecules from PLGA particles. It 
was also shown which aspects of the method of production 
coupled with molecular description of the formulation have 
the most impact on the release rate. Heuristic modeling tech- 
niques showed their usefulness in the complex and still not 
fully understood process of dissolution from PLGA particles. 
Moreover, they led to the mathematical formula (Equation 3) 
describing the process with predictive efficiency comparable 
to the ANN models. 

With predictive modeling, knowledge extraction was 
performed simultaneously as a part of the modeling pro- 
tocol. The major procedure employed feature selection 
techniques for reduction of input vector. A successful 
reduction from an initial 300 to 17 and eleven provided both 
information about the analyzed problem and appropriate 
generalization ability of MON-MLP and mathematical 
model, avoiding the well-known problem of "curse of 
dimensionality." Analysis of crucial variables confirmed 
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the dominating role of PLGA formulation characteristics; 
however, geometric characteristics of protein molecules and 
plasticizer ionic descriptor were also found to be impor- 
tant features. Analysis of the selected variables is a direct 
source of knowledge and a possible hint for creating a more 
general theory of how macromolecules are released from 
the PLGA particles. 

This work was possible thanks to the application of two 
key elements: 

• Marvin, a cheminformatic tool, which allowed direct 
computation of macromolecule descriptors - an ability 
not present in every cheminformatic package; 18 and 

• heuristic modeling tools for feature selection and predic- 
tive modeling integrating various types of information 
into the single model. 

The resulting models are versatile in terms of prediction of 
the unknown macromolecule release profiles from the PLGA 
particles, thus it might be useful for the future development of 
PLGA microparticles as a decision support system. 
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Table S3 The results of the formulation-to-formulation 



relRMSE (%) 



Formulation 


relRMSE (56) 




no 


MON-MLP 


GP model 




model 


(Equation 2) 


1 


9.56 


14.43 


2 


10.29 


1 1.70 


3 


22.14 


18.00 


4 


2.56 


1 1.74 


5 


22.63 


14.75 


6 


12.01 


19.47 


7 


10.72 


13.99 


8 


10.12 


13.87 


9 


1 1.42 


12.84 


10 


20.92 


21.92 


1 1 


17.73 


5.25 


12 


7.48 


8.78 


13 


17.47 


8.94 


14 


12.79 


23.37 


15 


3.12 


5.35 


16 


1.75 


9.75 


17 


16.10 


6.76 


18 


13.47 


8.76 


19 


1.83 


5.27 


20 


4.31 


8.38 


21 


12.97 


15.05 


22 


3.87 


7.57 


23 


6.26 


12.34 


24 


6.55 


6.99 


25 


5.35 


12.21 


26 


25.05 


13.12 


27 


14.48 


8.83 


28 


17.05 


14.67 


29 


1 1.20 


9.26 


30 


15.05 


7.01 


31 


36.03 


27.25 


32 


34.23 


26.38 


33 


8.21 


10.49 


34 


26.99 


25.33 
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Table S3 (Continued) 



Formulation 


reirxnot y/o) 




no 


MON-MLP 


GP model 




model 


(Equation 2) 


35 


6.85 


5.46 


36 


26.82 


22.38 


37 


28.55 


22.87 


38 


19.27 


24.45 


39 


8.78 


1 1.05 


40 


9.20 


9.04 


41 


18.09 


17.55 


42 


16.36 


15.15 


43 


10.71 


5.39 


44 


7.15 


10.71 


45 


20.51 


7.39 


46 


8.75 


32.83 


47 


15.23 


23.02 


48 


21.05 


20.61 


49 


4.70 


17.83 


50 


1 1.18 


6.56 


51 


10.90 


16.59 


52 


1 1.20 


18.85 


53 


10.85 


19.87 


54 


7.63 


10.20 


55 


25.31 


10.20 


57 


19.16 


24.17 


58 


24.05 


21.14 


59 


31.39 


29.31 


60 


28.53 


22.90 


61 


30.81 


5.76 


62 


18.73 


6.79 


63 


7.59 


8.59 


64 


13.37 


6.16 


65 


8.93 


7.73 


66 


4.61 


10.34 


67 


1 1.72 


4.76 


68 


14.25 


5.85 


Min 


1.75 


4.76 


Max 


36.03 


32.83 


Mean 


14.39 


13.75 
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